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ABSTRACT 

In the past few decades detailed observations of radio and X-rays emission from 
massive binary systems revealed a whole new physics present in such systems. Both 
thermal and non-thermal components of this emission indicate that most of the radi- 
ation at these bands originates in shocks. OB and WR stars present supersonic and 
massive winds that, when colliding, emit largely due to the free-free radiation. The 
non-thermal radio and X-ray emissions are due to synchrotron and inverse compton 
processes, respectively. In this case, magnetic fields are expected to play an important 
role on the emission distribution. In the past few years the modeling of the free-free 
and synchrotron emissions from massive binary systems have been based on purely 
hydrodynamical simulations, and ad hoc assumptions regarding the distribution of 
magnetic energy and the field geometry. In this work we provide the first full MHD 
numerical simulations of wind-wind collision in massive binary systems. We study the 
free-free emission characterizing its dependence on the stellar and orbital parameters. 
We also study self-consistently the evolution of the magnetic field at the shock region, 
obtaining also the synchrotron energy distribution integrated along different lines of 
sight. We show that the magnetic field in the shocks is larger than that obtained when 
the proportionality between B and the plasma density is assumed. Also, we show 
that the role of the synchrotron emission relative to the total radio emission has been 
underestimated. 

Key words: binaries: general - stars: winds - methods: numerical 



1 INTRODUCTION 

Massive stars (O/WR) are responsible for the emission of 
a large number of high energy photons that ionize most 
of the H and He present in the stellar envelope. Typically, 
the electron temperature can be of T ~ 10 4 — 10 5 K. The 
high density of the ionized envelope of WR stars results in 
large optically thick free-free emission (Falceta-Goncalves, 
Jatenco-Pereira & Abraham 2005, Abraham et al. 2005a,b). 
Wright & Barlow (1975) showed that a spherically symmet- 
ric expanding envelope, with a mass density profile p oc r~ 2 , 
results in a power law thermal radio emission with a positive 
spectral index a ~ 0.6, which is confirmed observationally 
for most WR stars (e.g. Williams 1996). However, some ob- 
jects, like sources with radio fluxes that peak at frequencies 
corresponding to T > 10 6 K, present significant smaller, or 
even negative, slopes (a < 0) for the power spectrum at 
these wavelengths (de Becker 2007). This emission is gen- 
erally attributed to synchrotron radiation of optically thin 
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plasma. The puzzle was to explain the physical conditions 
for a high temperature, but low density, plasma in the WR 
stellar envelope. 

Compared to individual WR stars, binaries present 
larger X-rays fluxes due to wind-wind shocks. The main 
source of this radiation is the radiative cooling, mostly free- 
free emission, from a hot plasma at T > 10 5 K. Theoretically, 
it is predicted that the interaction of the two supersonic 
winds of a WR+O binary system results in strong shocks. 
This scenario has been shown, both analytically and nu- 
merically, by several authors (e.g. Usov 1991, Stevens 1995, 
Walder 1998, Antokhin, Owocki & Brown 2004, Parkin & 
Pittard 2008, Pittard & Parkin 2010, and others) to be effi- 
cient in developing dense and hot plasma that can account 
for the observed X-ray emission. 

While free-free emission is produced by Coulombian col- 
lisions of charged particles, the non-thermal contribution is 
related to the relativistic electrons accelerated by magnetic 
fields. Eichler & Usov (1993) showed that, if magnetized, 
the strong shocks occurring at wind-wind interaction zones 
are responsible for the acceleration of charged particles to 
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relativistic speeds. Electrons may reach relativistic speeds 
in strong shocks due to the first order Fermi acceleration 
process. As the magnetic pressure also increases at shocks 
the synchrotron emission may dominate the total emission 
at radio wavelengths. 

Besides of providing a general picture, analytical models 
are unsuccessful in predicting the actual contribution of non- 
thermal sources because of the non-linear dynamical evolu- 
tion of the post-shock plasma. Under the very specific con- 
ditions of long period orbits, shock symmetry/homogeneity 
and slow cooling, analytical models are good approxima- 
tions. In a different scenario, where strong cooling is present 
or multidimensional description of variables is needed, nu- 
merical simulations should be used instead. Currently, nu- 
merical simulations represent the best method to fully un- 
derstand and describe the physics of complex wind- wind col- 
lision. 

Dougherty et al (2003) and Pittard et al. (2006) mod- 
elled the radio emission maps from massive binary systems 
based on the output of hydrodynamic numerical simulations. 
Pittard & Dougherty (2006) have also compared these mod- 
els to the observed features of WR140. Their simulations did 
not include the evolution of the magnetic fields. For that 
they assumed energy equipartition to obtain the distribu- 
tion of the magnetic energy density at the shock region. It is 
known however that the correlations between magnetic pres- 
sure, density and thermal energy are not linear (Burkhart 
et al 2009). Full MHD simulations are therefore mandatory. 

In this work we calculate a number of MHD numerical 
simulations of strong wind-wind collisions in order to study 
the distribution of the magnetic fields within the shock re- 
gion. Based on the self-consistent thermal and magnetic en- 
ergy distributions we also calculate the synthetic emission 
maps of synchrotron radiation from the simulated data, and 
compare them with those of thermal emission. In the fol- 
lowing section we review the basic physics of the wind-wind 
collision model and describe the numerical setup of the prob- 
lem. In Section 3, we present the main results, describe the 
gas and magnetic energy distributions at the shock region 
and obtain the synthetic maps for the radio emission. Fi- 
nally, the Conclusions are shown at Section 5. 



2 THE MODEL 

2.1 The physics of the wind-wind shocks 

In a massive binary system the two stellar winds are ex- 
pected to present high mass loss rates and supersonic veloc- 
ities. The interaction of these winds is a shock region formed 
between the two stars. The shock results in an increase of 
gas density and temperature, consequently also in emission 
at certain wavelengths (e.g. X-rays), being the total emission 
strongly dependent on the post-shock values of density and 
temperature. For the non-thermal emission, the magnetic 
field intensity and the population of relativistic particles is 
also needed. 

The shock front is formed at a distance n from the 
primary star given by: 

n = -^ (1) 

l + v 2 

where D is the distance between the two stars and r\ = 



M s Vs/MpVp. M p and M s are the mass loss rates of the pri- 
mary and the companion star, respectively, and v p and v s 
their asymptotic wind velocities. 

If magnetic fields are neglected, and the shock consid- 
ered adiabatic, the Rankine-Hugoniot (RH) jump conditions 
for density and temperature at pre (1) and post-shock (2) re- 
gions may be obtained by the basic hydrodynamic equations 
for mass continuity, fluid momentum and energy conserva- 
tion: 



PlVl — P2V2, 
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where p is the gas mass density, v the relative wind velocity, 
P the gas pressure and 7 the adiabatic exponent; indexes 1 
and 2 indicate pre (upstream) and post-shock (downstream) 
parameters. From Eqs.(2) - (4) we can calculate the jump 
relations for density and temperature: 
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where M = Vi ( / ykTi/pm H )~' s is the Mach number, p the 
molecular weight and mu the mass of the H atom. 

It is known however that this scenario is dramatically 
changed if cooling processes are considered. Cooling reduces 
the downwind gas pressure, shrinking the shocked gas into a 
thin layer. The plasma flow over this thin layer result in the 
transverse acceleration instability (Dgani, Walder & Nuss- 
baumer 1993, Dgani, Van Buren & Noriega-Crespo 1996) 
and/or the thin-layer instability (Vishniac 1994, Folini & 
Walder 2000). 

The transverse acceleration instability arises in ex- 
tremely thin shocks, which is a reasonable approximation 
for shocks with infinite Mach number and very strong cool- 
ing. As described in Dgani, Walder & Nussbaumer (1993), 
the ram pressure over the shocked fluid element depends 
on its distance to the line that connects the two stars. The 
instability arises as the perturbations in the tangential ve- 
locity along the shock surface result in a net torque that 
increases even more the angle of this disturbance (velocity 
vector) . 

On the other hand, the thin-layer instabilities, which 
lead to the growth of turbulent motions, may arise in strong 
shocks with fast cooling. The thin layer instability is under- 
stood as the post shock gas flow transports momentum away 
of the shock region. The slabs then produce substructures 
on scales of the shock thickness, which travel at basically 
the sound speed. More generally, the thin layer instability 
will grow at scales L < I < c s t, being L the shock surface 
thickness. The large broadening observed in spectral lines 
of certain binary systems indicates that turbulent motions 
may be present at the post-shock region (Falceta-Gongalves, 
Abraham & Jatenco-Pereira 2006). As pointed by Pittard 
(2007), clumpy pre-shock winds naturally leads to a turbu- 
lent post-shocked medium though the turbulent power, as 
well as the length scales of the largest eddies may vary sub- 
stantially from those excited by instabilities. 
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The instability growth rate, however, is reduced if the 
magnetic field is strong at the shock, because the magnetic 
pressure increases the shock width. Even in the strongly 
magnetized case, O-WR stars push the magnetic fields out- 
wards with their strong winds. At the shock, the magnetic 
fields also play a role in reducing the compression degree of 
the plasma, as described as follows. 

In magnetized shocks the momentum and energy jump 
conditions (Eqs.[3] and [4]) are changed, and another equa- 
tion, to account for the magnetic field evolution, is required. 
The magnetized RH jump conditions are given by: 

(7) 
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(10) 



where B is the intensity of the magnetic field component 
parallel to the shock front. 

For a weak shock, i.e. P2/P1 ~ 1, the above equations 
reduce to vi ~ 7c 2 + c? A , being c s = {kT/rriH) 1 the sound 
speed and c\ = B /(A-Kp) 1 ^ 2 the Alfven speed. Therefore, the 
magnetic field reduces the "strength" of the shock since part 
of the upstream kinetic energy is converted into magnetic 
despite of thermal energy. 

Notice that the increase in magnetic energy density is 
not related to any type of dynamo process here. Not in the 
adiabatic case at least. The magnetic field lines, frozen to 
the plasma, are trapped within the shock region by the com- 
pressed gas. For a shock with both sonic and alfvenic Macqj 
numbers M 3 > Ma 3> 1, the post-shock parallel component 
of the magnetic field is Am a ~ 4Sii l5 if magnetic diffusion is 
neglected, where B\\% is the intensity of the upstream mag- 
netic field component parallel to the shock surface. 

In the case of strong cooling the upstream thermal pres- 
sure is rapidly reduced, as energy is removed by radiation, 
and the shocked plasma is compressed even more in order to 
reach pressure equilibrium. In this situation both the density 
and the intensity of the parallel component of the magnetic 
field will increase more than a factor of ~ 4. 

At very high temperatures, such as in the wind-wind 
collisions, free-free emission is considered the main cooling 
mechanism. The free-free emissivity is given by 

P„ ff ~ 6.81x 10" 38 Z 2 n e mg s T~ i e~& (erg s~ 1 cm" 3 Hz" 1 ), (11) 

where gg is the averaged gaunt factor for free- free emission, 
and the free-free absorption coefficient 

nl ~ 3.7 x 10 8 Z 2 n e mgffT~K~ 3 (l - e"*r)(cm _1 ), (12) 

and, the optical depth for the free- free absorption is obtained 
by 



Tff 



K,,as. 



(13) 



1 the sonic and alfvenic Mach numbers are defined as the ratio of 
the wind velocity by the sound and Alfven speeds, respectively, 
i.e. M s = v/cs and Ma = v/ca, where ca = S(47rp) ' . 



2.2 Particle acceleration and the synchrotron 
emission 

As pointed by Eichler & Usov (1993), the conditions for ef- 
ficient particle acceleration at the shocks are: i) the shock 
being collisionless, i.e. the timescales for particle collisions 
must be large compared to ion cyclotron period, ii) the 
Ohmic damping rate of the Alfven modes to be small, and 
iii) the shock speeds must exceed the phase velocity of the 
whistler mode. These conditions are satisfied for: 
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UO'kms- 1 / Vl0 10 cm- 3 

Once the two conditions above are satisfied electrons 
will suffer multiple reflections and leave the region with in- 
creased energy. This process is balanced by the processes of 
energy loss, such as collisions with other ions, emission of 
radiation and the inverse Compton (IC) scattering. Typi- 
cally, synchrotron emission and IC scattering are the most 
important processes of energy loss for relativistic particles. 
The ratio between the IC and synchrotron losses is given by: 

7 - 77- -f*KN, (,-Lo) 

-ksyn UB 

where L/photons is the energy density of photons that will 
be scattered by the relativistic electrons, Ub = B 2 /8n is 
the magnetic energy density, and Pkn is the Klein-Nishina 
correction factor. 

For the case in which the IC is dominant, i.e. U p h > 
Ub, the maximum energy of electrons accelerated at the 
magnetized shock is -E ma x = 7max?noc 2 , being: 
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where tl is the Larmor radius, A m f p the mean free path, and 
o"t the cross-section for Thompson scattering. 

For typical parameters of wind-wind collisions in mas- 
sive binary systems it is possible to obtain particles with 
very high energies (7 > 10 2 ). The back reaction of particles 
in the shock structure is of course important and may re- 
duce the average acceleration (Pittard & Dougherty 2006). 
In this case, the number of particles accelerated to very high 
energies will be reduced. In any case, a small fraction (x) of 
the total electrons will be accelerated to relativistic speeds, 
modifying the tail of the Maxwellian distribution of thermal 
particles. In the linear regime, the population of relativistic 
electrons is given by the power law (Bell 1978): 



/(p) ~ c P 



_ C + 2 



(18) 



where C is the normalization constant, p the particle linear 
momentum normalized by m e c and £ is the compression 
index of the shock p/po- Typically, f ~ 4 for strong adiabatic 
shocks, resulting in a power-law index s = (f +2)(f — 1) ~ 2. 
The stochastic acceleration occurs efficiently for elec- 
trons with speeds larger than the sound speed. The frac- 
tion of relativistic to the total number of electrons then de- 
pends on the minimum momentum necessary for the parti- 
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Figure 1. Logarithm of the plasma density (top) and temperature (bottom) obtained for the adiabatic case (left) and the model with 
radiative cooling (right). 
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Figure 2. Logarithm of the magnetic field intensity for the adiabatic case (left) and the model with radiative cooling (right). The scale 
is normalized by the magnetic field intensity at the surface of the stars. 



cles to be accelerated (p™j n ), and may be obtained from the 
Maxwellian distribution (<^>(p)) as follows: 



4>(p)dp ■ 
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which gives X ~ 10~ 4 - 10~ 3 , for p™j n ~ 3.3-3.6. When the 
first order Fermi acceleration process is dominant, and low 
collisionality of particles is assumed, % may be even larger 
reaching few percent (Ellison & Eichler 1985). In this work, 
in the following sections, we use x = 10~ 3 . 
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Figure 3. Emission maps of free-free (top) and synchrotron (bottom) emission at v = 1GHz calculated for the adiabatic (left) and 
cooled models (right). The scale is normalized by the maximum emission obtained for the free- free emission from the adiabatic case. 



The relativistic electrons will then emit synchrotron ra- 
diation. For a power-law distribution of electron momenta 
the synchrotron emission is given by: 

(erg s _1 cm _3 Hz _1 )(20) 



S— 1 nl+O, 



with a — (s — l)/2 being the power-law index of the syn- 
chrotron emission distribution. 

The emission at lower frequencies is quenched by the 
Tsytovitch-Razin effect, i.e. within a plasma with refractive 
index lower than unity the phase velocity of electromagnetic 
waves is larger than c. As a result, the effective Lorentz fac- 
tor is largely reduced for frequencies in which the refractive 
index is much lower than unity. Following Dougherty et al. 
(2003), a simple way of modeling the Tsytovitch-Razin effect 
in the calculation of the synchrotron emission is assuming 
an exponential decay below the cut-off frequency, 

psyn,Razin = ps y n e -^ ) (21) 

where vr ~ v v Jvb-, being v v and vb the plasma and cy- 
clotron frequencies respectively. 

The radial and azimuthal profiles of gas density and 
temperature at the shock region are not simple, therefore the 
proper analytical calculation of total emissions for free-free 
and synchrotron radiation is very complex. Also, the evolu- 
tion of the gas for non-adiabatic shocks (i.e. when cooling is 
fast is a non-linear process and its description by analytical 



approximations is not possible. Numerical MHD simulations 
are mandatory for the proper determination of the density, 
temperature and magnetic field intensity profiles. 



2.3 Numerical Simulations 

The simulations were performed solving the set of ideal 
MHD equations, in conservative form, as follows: 

(22) 
= f , (23) 

(24) 

(25) 

where p, v and p are the plasma density, velocity and 
pressure, respectively, B is the magnetic field and f rep- 
resents the external source terms (e.g. gravity of both 
stars) is the magnetic field. The equations are solved using 
a second-order-accurate and non-oscillatory scheme, with 
open boundaries. 

The set of equations is only complete with an equation 
of state, or with an explicit equation for the evolution of the 
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energy. In this work we studied two different cases, one is 
assuming that the plasma is adiabatic and in the second the 
radiative cooling is included in the calculations. 

In the adiabatic case, the set of equations is closed with, 



Pocp 1 



(26) 



being 7 = 5/3 the adiabatic constant. 

For the radiative cooling the set of equations is closed 
calculating 

riP 1 

7» = — n2A(n (27) 

after each timestep, where n is the number density, P is the 
gas pressure, 7 the adiabatic constant and A(T) is the inter- 
polation function from an electron cooling efficiency table 
for an optically thin gas (Gnat & Sternberg 2007). 

The initial setup for the stellar winds was chosen based 
on the interesting case of 77 Carinae. The primary star mass 
loss rate is assumed as M\ ~ 2 x 10 -4 M Q yr _1 with a 
stellar wind terminal velocity vi ~ 700 km s~ . For the 
secondary star we used M2 — 10 M© yr -1 and Vi — 3 x 
10 3 km s , respectively (Pittard & Corcoran 2002, Falceta- 
Gongalves, Jatenco-Pereira & Abraham 2005, Abraham & 
Falceta-Gongalves 2007). The distance between the stars is 
set = 0.8 AU, which corresponds to the periastron passage of 
an orbit with a major semi-axis a — 15AU and eccentricity 
e = 0.95. The simulated box is set as a fixed eulerian grid of 
512 3 cells, being L — 5 x 10 15 cm its size in each direction. 
For the sake of simplicity we used solar abundances for the 
chemical composition of the gao 

The magnetic field is initially set as a dipole at the 
surface of both stars with polar intensity Bo = 1G, cor- 
responding to a thermal to magnetic pressure ratio fit = 
P thermal /Pmag ~ 20 for the primary star, and 02 ~ 5 for the 
secondary. The orientation of the magnetic dipoles in the ini- 
tial setup is irrelevant because of the high f3 values assumed 
in the calculations. For the chosen initial magnetic field in- 
tensities the wind kinetic and thermal pressures dominate 
the dynamics of the plasma. The magnetic field plays minor 
role and is dragged by the wind flows, to an approximately 
radial distribution (B oc r~ 2 ), as it expands. 



3 RESULTS 

As mentioned above, we run two different models, one con- 
sidering the adiabatic case and another taking into account 
the radiative cooling in the evolution of the shocked gas. 
The simulations were carried out until t — 2.5 x 10 6 s. The 
profiles of density, temperature and magnetic field of the fi- 
nal data cube are presented below. In Figure 1 we show the 
central slices of density (up) and temperature (bottom) for 
each model, the adiabatic (left) and with radiative cooling 
(right). 

2 The chemical composition of the winds is relevant when calcu- 
lating the degree of ionization, fraction of electrons to ions and 
average charge of ions, for instance. Eta Carinae (LBV star) and 
its companion, possibly a WR star, present different abundances. 
However, the code used in these calculations is single-fluid in the 
sense that we cannot separate the contributions from each wind 
after partial mixing at the shock region. 



The adiabatic case presents the standard prominent 
features of wind-wind collisions, such as the primary and 
secondary shock surfaces, and the contact discontinuity be- 
tween them. As expected, the increase in density at the post- 
shock of each surface is ~ 4. The primary star wind is denser 
than the one of the secondary star. Therefore, the densest 
emitting plasma is located between the primary shock sur- 
face and the contact discontinuity, being the density almost 
two orders of magnitude lower at the secondary shock layer. 
On the other hand, the temperature jump at the shock is 
proportional to the Mach number and is larger at the sec- 
ondary shock layer. The temperature peaks ~ 5 x 10 7 K at 
the secondary shock surface while it reaches ~ 10 6 at the 
primary. 

When the radiative losses are taken into account this 
picture changes dramatically. The energy loss occurs rapidly, 
in a timescale shorter than the expansion of the heated gas 
along the shock surfaces, i.e. 



T s] 
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layer 



nA(T shock ) v, 



(28) 



being the cooling time t coo i ~ T/nA(T), the dynamical 
time Tdy n ~ l/v, /layer the typical length scale of the shock, 
and «ex P the flow speed at the downstream. 

In such a case the total pressure within the shock sur- 
faces is largely reduced, and the external pressure of the 
winds causes the contraction of the region. The shock re- 
gion, thinner compared to the adiabatic case, is then sub- 
ject to the thin-layer instability, clearly seen at both the 
density and temperature maps. The most prominent differ- 
ence here is the absence of the three discontinuities seen 
in the adiabatic case, since the shocked gas is now mixed, 
compressed in a smaller region and turbulent. Due to the 
compression, the density is also larger compared to the adi- 
abatic case (n max ~ 10 12 cm -3 , and the jump of density in 
the shock 3> 4.), as predicted in the analytical approach 
by Falceta-Gongalves et al. (2005). Due to the cooling af- 
ter the shock, the temperature of the downstream plasma 
is smaller compared to the adiabatic case, reaching maxima 
around ~ 10 6 K. 

The magnetic field intensity at the shocks, for both 
models, is shown in Figure 2. As mentioned above, the mag- 
netic field lines are basically radial at the flows up to the 
shocks. The flow however is not perpendicular to the shock 
(except for the line intercepting the two stars) and there is 
always a field component that will be amplified. In the adi- 
abatic case the amplification of the magnetic field intensity 
is larger at the contact discontinuity. At the shock surfaces 
the increase in the field intensity (the magnetic field com- 
ponent parallel to the shock surface only) is around ~ 4, 
as expected, since its jump condition is proportional to the 
gas density. However, the magnetic energy density increases 
towards the contact discontinuity. The reason for this is the 
orientation of the field lines. As the downstream plasma 
propagates towards the contact discontinuity it drags the 
magnetic field lines with it. When the material reaches the 
contact discontinuity starts to flow along it. The field lines 
also get aligned with the contact discontinuity. At this stage, 
since the flow is parallel to the field lines, the magnetic field 
lines are not dragged outwards. Therefore, two processes 
take place; one is field lines being dragged into the shock 
region and second the flow running outwards leaving the 
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field lines piled up at the contact discontinuity. This effect 
causes a gradual increase of the magnetic pressure up to the 
cquipartition. 

For the model with radiative losses the degree of com- 
pression is large and, as a consequence, the magnetic pres- 
sure is also increased. As more gas is confined in the shock 
region the field lines, frozen to the plasma, accumulate at 
the shock region. The magnetic field intensity in this model 
is approximately two orders of magnitude larger than in the 
adiabatic case (Bmax ~ 10 2 -Bmax)- Due to the turbulent na- 
ture of the post-shocked flow the morphology of the field 
lines is also complex, in contrast to the uniform distribution 
of the field lines in the adiabatic case. In this model, the 
magnetic energy is concentrated in filaments more or less 
well distributed along the shock region, despite of the con- 
centration of magnetic energy seen in the adiabatic model, 
at the shock apex. 

These differences in density, temperature and magnetic 
energy distributions for the adiabatic and non-adiabatic 
models play a major role in the spectral energy distribution 
of the shock emission. 
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Figure 5. Total emission calculated for the adiabatic case: free- 
free (solid) and synchrotron (dashed), and for the case where the 
cooling module was used: free-free (dotted) and synchrotron (dot- 
dashed). All fluxes are normalized by the total free- free emission 
at 1GHz of the adiabatic case. 



3.1 Synthetic radio maps 

Considering the amplifications seen in the magnetic pres- 
sure at the shock and its complex distribution along it, the 
contributions of free-free and synchrotron emissions to the 
total radio emission of close binary systems may change sig- 
nificantly. 

To test this, we calculate the emission maps of both free- 
free and synchrotron emissions in the range of 1 to 100GHz. 
In order to obtain the total emission we perform the calcu- 
lation of the free-free and synchrotron emissivities (Eqs. 11 
and 20) at each cell of the cube, as well as the absorption 
coefficients and the Tsytovitch-Razin effect in the reduction 
of the synchrotron emission. 

The radiative transfer is then calculated for a chosen 
line of sight. In Figure 3, we show the free-free and syn- 
chrotron emissions calculated for the adiabatic and non- 
adiabatic cases, along z-direction (i.e. perpendicular to the 
plane shown in Figs. 1 and 2). The integrated emissions are 
obtained as the summation below: 



Pu(i,j) 



V x 



E 
fc=i 



P v {i,j,k)e- T ^' k \ 



(29) 



being V the volume of each cell in real units, and the optical 
depth obtained as: 



T„(i,j,k) = V3 



k = k 



K v {l,J, k), 



(30) 



The left panels of Figure 3 represent the free-free and 
synchrotron radio emissions for the adiabatic model, calcu- 
lated at 1GHz. The color bars indicate the fluxes normalized 
by the maximum value obtained for the free- free emission. 
The synchrotron emission is lower compared to the thermal 
source; however the most interesting feature of these panels 
is the location of the flux maxima. The free-free emission 
follows more or less the density distribution, being maxi- 
mum at the shock apex, decreasing mildly along the primary 
star shock. Also, the distribution along the map is smooth. 



The opposite side of the shock region, i.e. the secondary star 
shock, is much less relevant for the thermal emission. On the 
other hand, for the non-thermal emission, most of the flux is 
originated in the contact discontinuity where the magnetic 
field is actually amplified. This result is in disagreement to 
the assumption made by e.g. Dougherty et al. (2003) that 
the magnetic energy density follows linearly the distribution 
of the thermal pressure. This fact severely affects the direct 
comparison of synthetic thermal and non-thermal emission 
maps to the spatially resolved observations. 

For the case with radiative losses taken into account 
(right panels in Figure 3), the scenario is completely dif- 
ferent. Again the fluxes are normalized by the maximum 
emission of the thermal emission. In this case, the peak of 
the non-thermal emission is about two orders of magnitude 
larger than the maximum thermal flux. The clumpy appear- 
ance of the synchrotron emission is visualized for different 
angles in Figure 4. 

The spectral energy distribution in the range of 1 — 100 
GHz is shown in Figure 5. Here we integrate the fluxes of all 
pixels of the synthetic emission maps obtained for each fre- 
quency. For the case of r^Car, the free-free emission is domi- 
nant for the adiabatic case, being the synchrotron emission 
several orders of magnitude lower. The optical depth of the 
latter is low in this case, and the spectral index for the non- 
thermal source is negative (a ~ —0.2). The break at lower 
frequencies is due to the Tsytovitch-Razin effect. 

For the non-adibatic model, the contribution of each 
mechanism is quite similar, with an inversion of impor- 
tance around 10GHz. For instance, considering the physi- 
cal properties at the shock apex, the gas number density is 
n ~ 10 1 2cm -3 and the magnetic field strength is of order 
of 2kG. In this particular case, the cutoff frequency for the 
Tsytovitch-Razin effect is of order of 14GHz. Obviously, it 
varies at each cell due to the changes in the local properties, 
being the cutoff of most of the shocked region around 1GHz. 

At low frequencies synchrotron emission is dominant 
(as seen in Figure 3), but the thermal radiation dominates 
at larger frequencies. The positive slope for the synchrotron 
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Figure 4. Visualizations of the synchtrotron emission maps for different orientations of the LOS. 



emission (a ~ 0.6 — 0.7) reveals the importance of absorp- 
tion in this case. Pittard (2010) presented an extensive study 
of the properties of thermal radio emission from O+O bi- 
nary systems in which he showed the increase in the thermal 
emission for radiative cooling shocks, as the free- free emis- 
sion relates to temperature and density as e// oc n 2 T ' 2 . As 
the shock cools, the larger density results in stronger emis- 
sion. However, in the case presented in this work absorption 
becomes more important. The absorption results in lower 
emission when compared to the adiabatic case. As the ab- 
sorption coefficient a// scales as v~ 3 , we expect a major ef- 
fect at lower frequencies, as we showed in Figure 3 where the 
maps were obtained for 1GHz. The difference in the spectral 
indices for both models reveal the transition from optically 
thin to optically thick dominated free-free emission at the 
range of 10 - 100 GHz. 

As we showed above, the switch between thermal- 
dominated to synchrotron-dominated emissions in binary 
systems strongly depends on the magnetic field intensity and 
density distributions along the shock. The studied system of 
??Car is a case of very close binary system (< 1AU at perias- 
tron) and of winds with large mass-loss rates. Short period 
or highly eccentric orbits are also found in many WR+O and 
O+O systems, so we expect to observe a similar behavior in 
the spectral energy distribution of these systems as well. For 
the cases of long period (not eccentric) orbits these conclu- 
sions may not stand (see Pittard 2010). Even though not 
simulating these specific conditions, we may conjecture that 
an order of magnitude lower mass-loss rate of the primary 
star, and an increase in one order of magnitude in the stel- 
lar magnetic field would result in a dominant non-thermal 
emission even for the adiabatic case. 



4 DISCUSSION 

In this work we presented the first full MHD simulations 
of a wind-wind collision. We chose the stellar and orbital 
parameters of the binary system of r;Car in order to study 
the fraction of the free-free and synchrotron emissions. We 
focused at the physical conditions and the radio emission, 
in the range of 1 — 100 GHz, at the periastron passage. 

Several previous works aimed the study of thermal and 
non-thermal emissions at these frequencies based on hydro- 
dynamical numerical simulations (e.g. Dougherty et al. 2003, 
Pittard & Dougherty 2005, Pittard et al. 2006). In these pa- 
pers, the magnetic field was assumed to behave in a linear 
way, being the magnetic energy density proportional to the 



thermal pressure everywhere in the cube. Important conclu- 
sions made, such as the importance of the synchrotron self- 
absorption and the Tsytovitch-Razin effect in the observed 
spectral brake, were strictly dependent on this assumption. 

The most important result in this work is the self- 
consistent determination of the evolution of the magnetic 
field in the shock, for adiabatic and non-adiabatic cases. 
In both cases, the magnetic pressure does not follow the 
thermal pressure linearly - as assumed in previous works. 
Regarding the special case of r/Car, during the periastron 
passage, the conclusions regarding spectral features diverge 
from the ones listed by other authors - at least in the range 
of frequencies studied in this work. 

The observed break in the spectra at 10GHz for the 
adiabatic case is due to the thermal absorption. The cut- 
off frequency for the Tsytovitch-Razin effect occurs for the 
simulated cubes approximately at vr ~ 1.8GHz. 

Pittard et al. (2006) showed the importance of Inverse 
Compton effect in the cooling of the central regions of the 
shock (applied to WR 147). From our simulations, we show 
it is not the case for the non-adiabatic case. In that model, 
the synchrotron losses are dominant 

LlC _ ^photons 



10" 



(31) 



where [/photons an d Ub represent the photons and magnetic 
energy density, respectively. For the adiabatic case the ratio 
above is about unity. Obviously, the magnetic field intensity 
at the stellar surfaces has been assumed in this work to be 
1G. In a different situation, with much smaller magnetic 
energies, the Inverse Compton effect would be dominant. 

Finally, we showed that the peaks of the synchrotron 
and thermal emissions are not located at the same region 
of the shock. For the adiabatic case, the synchrotron emis- 
sion occurs mostly close to the contact discontinuity, where 
the magnetic energy density is larger. The thermal free-free 
emission on the other hand has its peak in the primary wind 
shock surface, due to the increased plasma density in the 
region. Therefore the direct comparison of the morphology 
of synthetic emission maps with those obtained from spa- 
tially resolved radio observations is not straightforward (e.g. 
WR 146 and WR 147). This is not an issue if the radiative 
cooling is fast compared to the dynamical timescales. In this 
situation, the shock region shrinks into turbulent filaments 
and the thermal and non-thermal emissions are overlapped. 
Fast cooling is expected to occur in short-period binary sys- 
tems, or during the periastron passage of long-period sys- 
tems but with very eccentric orbits (e.g. ?yCar). 

It is worth mentioning that the current model is still 
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limited in terms of several interesting processes that may be 
important in such systems. For instance, we did not treat 
the orbital motion of the stars. It is well known that short 
period, or highly eccentric systems, present a much more 
complex shock structure. Also, clumpy stellar winds may 
substantially change the distribution of density and velocity 
fields. In this case, the magnetic fields may vary in both com- 
plexity and intensity, even in the adiabatic case. The com- 
plexity of the magnetic field lines is related to the magnetic 
reconnection rate. In the models presented in this work, we 
did not treat the resistivity explicitly, but numerical resistiv- 
ity is always operating and diffusion of field lines and small 
scale reconnection still occurs, though its effects in heating 
the flow and accelerating particles, on the other hand, are 
not present. All these effects will be subject of study in fu- 
ture works. 



4.1 Magnetic reconnection and particle 
acceleration 

It has been proposed that magnetic field reconnection is in- 
creased due to small scale turbulent motions (Kowal, Lazar- 
ian & Vishniac 2009). At the current sheets turbulent recon- 
nection events of particle acceleration take place (Kowal, de 
Gouveia Dal Pino & Lazarian 2011). 

From the models presented in this work we expect par- 
ticle to be accelerated to relativistic speeds, followed by 
strong radio synchrotron and gamma ray emissions, in the 
shocked regions of binary systems due to turbulent reconnec- 
tion, rather than first order Fermi acceleration. Also, in this 
sense, we expect systems with large eccentricities to present 
variability in both particle acceleration efficiency and non- 
thermal emissions. As the secondary star moves away from 
the primary towards apastron the wind density at the shock 
decrease, as well as the cooling. As a consequence, the com- 
pression of the magnetic field lines, growth of instabilities 
and turbulence, and the reconnection rates will all decrease. 
As the secondary is closer to the primary, near periastron, 
the opposite would occur. 



• the maxima of synchrotron emission occur along the 
region between the contact discontinuity and the secondary 
star shock surface. The free- free emission, on the other hand, 
occurs mostly between the contact discontinuity and the pri- 
mary star shock surface. 

• if the radiative cooling is fast, as expected in close bi- 
nary systems or during periastron passage of long period but 
eccentric systems, the emission maps of thermal and non- 
thermal emissions are clumpy, due to the turbulent nature 
of the flow, and are overlapped. 

• the free-free emission dominates in the adiabatic case, 
and at the high frequencies range in the non-adiabatic case. 
The synchrotron emission is larger at low frequencies (1 — 
10GHz). 

• effects of absorption are important for the non-adiabatic 
case, in which the synchrotron spectral slope is positive. 

• for most of the massive binary systems, with mass-loss 
rates lower than those in ??Car, the situation is expected to 
be different, being the synchrotron emission the most impor- 
tant in the range of radio frequencies studied in this work. 

With more computational resources, we plan to extend 
this work in the near future by providing more simulations 
with different physical properties in order to validate these 
conclusions, or define the range in which they are still valid. 
Also, we plan to implement a more consistent method, inte- 
grated to the numerical calculations of the MHD equations, 
to determine the properties of the relativistic particles, such 
as fraction to total and maximum energies, which are impor- 
tant for the proper estimation of the synchrotron emission 
in these systems. 



ACKNOWLEDGMENTS 

The authors thank the Brazilian agencies FAPESP (no. 
2011/12909-8) and CNPq for financial support. We also 
thank the anonymous referee for the comments that helped 
improving this paper. 



5 CONCLUSIONS 

In this work we provide the results of the first full MHD 
numerical simulations of wind-wind collisions in binary sys- 
tems. We considered adiabatic equation of state and radia- 
tive cooling function as two different closures for the system 
of MHD equations. We evolved the basic physical parame- 
ters of the plasma and calculated the synthetic observable 
maps of free-free and synchrotron emissions in the range of 
I — 100GHz. For this work we used the physical parameters 
of the massive binary system of ??Car during its periastron 
passage. The main conclusions of the analysis made are, for 
the set of initial parameters chosen in this work: 

• the magnetic field evolution at the shock is not linearly 
correlated to the density, or thermal pressure. In the adi- 
abatic case, the magnetic field is amplified mostly close to 
the contact discontinuity. In the non-adiabatic case the flow 
is turbulent due to the growth of the thin-layer instabilities. 
The magnetic field amplification is 2 — 3 orders of magni- 
tude larger and its distribution is chaotic as it is frozen to 
the turbulent flow. 
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